library(forecast)
library(fpp2)
## Loading required package: fma
## 
## Attaching package: 'fma'
## The following object is masked from 'package:plyr':
## 
##     ozone
## Loading required package: expsmooth
library(tidyverse)
library(tseries)
detach("package:dplyr", unload=TRUE)
## Warning: 'dplyr' namespace cannot be unloaded:
##   namespace 'dplyr' is imported by 'tigris', 'broom', 'tidycensus' so cannot be unloaded
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:xts':
## 
##     first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(scales)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
library(GGally)
## 
## Attaching package: 'GGally'
## The following object is masked from 'package:dplyr':
## 
##     nasa
## The following object is masked from 'package:fma':
## 
##     pigs
library(e1071)
library(fpp2)
theme_set(theme_minimal())

3.1 Some simple forecasting methods

Average method

\[\hat{y}_{T+h|T} = (y_{1} + \dots + y_{T})/T\]

Naive method

\[\hat{y}_{T+h|T} = y_{t}\]

Seasonal Naive method

\[\hat{y}_{T+h|T} = y_{T+h-m(k+1)}\]

For example, with monthly data, the forecast for all future February values is equal to the last observed February value.

Drifted method

\[\hat{y}_{T+h|T} = y_T + \frac{h}{T-1}\sum_{t=2}^T(y_t-y_{t-1}) = y_T +h \frac{y_T - y_1}{T-1}\] This is equivalent to drawing a line between the first and last observations, and extrapolating it into the future.

x <- ts(1:20, start = 1990, frequency = 4)

autoplot(x)

autoplot(naive(x, 5))

autoplot(rwf(x, 4))

autoplot(snaive(x, 4))

beer2 <- window(ausbeer, start = 1992, end = c(2007, 4))

autoplot(beer2) +
  autolayer(meanf(beer2, h = 11),
            series = "Mean", PI = F)+
  autolayer(naive(beer2, h = 11),
            series = "Naive", PI = F) +
  autolayer(snaive(beer2, h = 11),
            series = "Seasonal Naive", PI = F) +
  autolayer(rwf(beer2, h = 11, drift = T),
            series = "Naive with Drift", PI = F) +
  labs(x = "Year", y = "Megalitres",
       title = "Foecasts for quarterly beer production")

These simple method will be consider as benchmarks. all new methods will be compare to these mmethod, to see whether method are better or not.

3.2 Transformation and adjustments

Calender adjustments

autoplot(milk)

dframe <- cbind(monthly = milk,
                DailyAverage = milk / monthdays(milk))

autoplot(dframe, facet = T) +
  labs(x = "year", y = "Pounds",
       title = "Milk production per cow")

Population adjustments

Inflation Adjustments

Mathematical Transformations

log transformations, power transformations, Box Cox transformations

lambda <- BoxCox.lambda(elec)
lambda
## [1] 0.2654076
autoplot(BoxCox(elec,lambda))

kpss.test(elec)
## Warning in kpss.test(elec): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  elec
## KPSS Level = 7.9424, Truncation lag parameter = 5, p-value = 0.01
dframe <- cbind(month = elec,
                DailyAverage = elec / monthdays(elec))

autoplot(dframe, facets = T) +
  labs(x = "year")

Bias adjustments

fc <- rwf(eggs, drift = T, lambda = 0, h = 50, level = 80)
fc2 <- rwf(eggs, drift = T, lambda = 0, h = 50, level = 80, biasadj = T)

autoplot(eggs) +
  autolayer(fc, series = "Simple back transformation") +
  autolayer(fc2, series = "Bias adjusted", PI = F) +
  guides(color = guide_legend(title = "Forecast"))

3.3 Residual Diagnostics

Check residual normality like linear regression does 1. residual uncorrelated 2. Residuals has a zero mean 3. Residuals have constant variance 4. The residuals are normally distributed.

Example: Forecasting the Google dialy Closing stock price.

\[e_{t} = y_t - \hat{y_t} = y_t - y_{t-1}\]

autoplot(goog200) +
  labs(x = "Day", y = "Close Price (US$)",
       title = "Google Stock (daily ending 6 December 2013)")

we can chekc the residual obtain from forecasting this series using the naive method

# plot the naive forecasting

goog_naive <- naive(goog200, 100)

autoplot(goog200) +
  autolayer(goog_naive, se = F, series = "Naive") +
  labs(x = "Day", y = "Close Price (US$)",
       title = "Google Stock (daily ending 6 December 2013)")

plot the residuals

res <- residuals(goog_naive)

t.test(res)
## 
##  One Sample t-test
## 
## data:  res
## t = 1.5892, df = 198, p-value = 0.1136
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1678207  1.5612704
## sample estimates:
## mean of x 
## 0.6967249
autoplot(res) +
  labs(x = "Day", y = "",
       title = "Residuals from naive method")

gghistogram(res) +
  labs(title = "Histogram of residuals") +
  stat_function(fun = dnorm,
                color = "red",
                args = list(mean = mean(res, na.rm = T),
                            sd = sd(res, na.rm = T)))
## Warning: Removed 1 rows containing non-finite values (stat_bin).

ggAcf(res) +
  labs(title = "ACF of residuals")

These graphs show that the naive method produces forecasts that appear to account for all avaiable info. The mean of the residuals is calose to zero, in t test result. There is no correlation in the residual series, in acf test no value exceed \[\pm2/\sqrt{n}\], n is the length of the Time Series. The times plot of the residuals show that the variation of teh residuals stay much the same across the historical data, aprat from the outlier, and therefore residual can be treated as constant.

Consequently, forecasts from naive method is a good fit, but the prediction interval that are computed assuming a normal distributuon may be inaccurate.

Portmanteau tests for autocorrelation

# lag = h and fitdf = K

Box.test(res, lag = 10, fitdf = 0)
## 
##  Box-Pierce test
## 
## data:  res
## X-squared = 10.611, df = 10, p-value = 0.3886
Box.test(res, lag = 10, fitdf = 0, type = "Lj")
## 
##  Box-Ljung test
## 
## data:  res
## X-squared = 11.031, df = 10, p-value = 0.3551

Since both Q and \(Q^*\), the result are no significatnt. Thus, we can clude that the residuals are not distriguishable from a white noise sieres.

All those methods for checking residuals are conveniently in R function checkresiduals()

checkresiduals(naive(goog200))

## 
##  Ljung-Box test
## 
## data:  Residuals from Naive method
## Q* = 11.031, df = 10, p-value = 0.3551
## 
## Model df: 0.   Total lags used: 10

3.4 Evaluating forecast accuracy

Training and test sets

start(ausbeer)
## [1] 1956    1
end(ausbeer)
## [1] 2010    2
window(ausbeer, start = 1995)
##      Qtr1 Qtr2 Qtr3 Qtr4
## 1995  426  408  416  520
## 1996  409  398  398  507
## 1997  432  398  406  526
## 1998  428  397  403  517
## 1999  435  383  424  521
## 2000  421  402  414  500
## 2001  451  380  416  492
## 2002  428  408  406  506
## 2003  435  380  421  490
## 2004  435  390  412  454
## 2005  416  403  408  482
## 2006  438  386  405  491
## 2007  427  383  394  473
## 2008  420  390  410  488
## 2009  415  398  419  488
## 2010  414  374
window(ausbeer, start = length(ausbeer) - 4*5)
## Warning in window.default(x, ...): 'start' value not changed
##      Qtr1 Qtr2 Qtr3 Qtr4
## 1956  284  213  227  308
## 1957  262  228  236  320
## 1958  272  233  237  313
## 1959  261  227  250  314
## 1960  286  227  260  311
## 1961  295  233  257  339
## 1962  279  250  270  346
## 1963  294  255  278  363
## 1964  313  273  300  370
## 1965  331  288  306  386
## 1966  335  288  308  402
## 1967  353  316  325  405
## 1968  393  319  327  442
## 1969  383  332  361  446
## 1970  387  357  374  466
## 1971  410  370  379  487
## 1972  419  378  393  506
## 1973  458  387  427  565
## 1974  465  445  450  556
## 1975  500  452  435  554
## 1976  510  433  453  548
## 1977  486  453  457  566
## 1978  515  464  431  588
## 1979  503  443  448  555
## 1980  513  427  473  526
## 1981  548  440  469  575
## 1982  493  433  480  576
## 1983  475  405  435  535
## 1984  453  430  417  552
## 1985  464  417  423  554
## 1986  459  428  429  534
## 1987  481  416  440  538
## 1988  474  440  447  598
## 1989  467  439  446  567
## 1990  485  441  429  599
## 1991  464  424  436  574
## 1992  443  410  420  532
## 1993  433  421  410  512
## 1994  449  381  423  531
## 1995  426  408  416  520
## 1996  409  398  398  507
## 1997  432  398  406  526
## 1998  428  397  403  517
## 1999  435  383  424  521
## 2000  421  402  414  500
## 2001  451  380  416  492
## 2002  428  408  406  506
## 2003  435  380  421  490
## 2004  435  390  412  454
## 2005  416  403  408  482
## 2006  438  386  405  491
## 2007  427  383  394  473
## 2008  420  390  410  488
## 2009  415  398  419  488
## 2010  414  374
tail(ausbeer, 4*5)
##      Qtr1 Qtr2 Qtr3 Qtr4
## 2005            408  482
## 2006  438  386  405  491
## 2007  427  383  394  473
## 2008  420  390  410  488
## 2009  415  398  419  488
## 2010  414  374

Forecast errors

Forecast error is the difference between an observation and its foecast

\[e_{T+h} = y_{T+h} - \hat{y}_{T+h|T}\] forecast error different from residuals in two ways

  1. Residuals are calcuated on the training set while forecast errors are calcuated on the test set.
  2. Residuals are based on on-step forecast while forecast eoors can involve multi-step forecasts.

we can measure forecast accuracy by summarising the forecast erro in different ways.

Scale-dependent erros

Mean absolute error: MAE = mean\(|e_t|\) Root mean squared error: RMSE = \(\sqrt{mean(e_t^2)}\)

When comparing forecast methods applied to a single time series, or to serveral time series with the same units, the MAE is popular as it it easy to both understand and compute. A forecast method that minise the MAE will lead to forecasts of the median, while minimising the RMSE will lead to forecast of the mean. Consequently, the RMSE is also widely used, depsite being more difficult to interpret.

Percentage errors

\(p_t = 100e_t/y_t\), mean unit-free

Mean absolute percentage error: MAPE = mean(\(|p_t|\)) Potential drawdown, being infinite or undefined if \(y_t = 0\) for any t in the period of interest.

Scaled errors

Scaled errors as an alternative to using percentage errors when comparing forecast accuracy across series with different units.

Examples

Seaonal dataset

beer2 <- window(ausbeer,start=1992,end=c(2007,4))

beerfit1 <- meanf(beer2,h=10)
beerfit2 <- rwf(beer2,h=10)
beerfit3 <- snaive(beer2,h=10)

autoplot(window(ausbeer, start=1992)) +
  autolayer(beerfit1, series="Mean", PI=FALSE) +
  autolayer(beerfit2, series="Naïve", PI=FALSE) +
  autolayer(beerfit3, series="Seasonal naïve", PI=FALSE) +
  labs(x = "year", y = "Megalitres", title = "Forecasts for quarterly beer production")

beer3 <- window(ausbeer, start = 2008)

accuracy(beerfit1, beer3)
##                   ME     RMSE      MAE        MPE     MAPE     MASE
## Training set   0.000 43.62858 35.23438 -0.9365102 7.886776 2.463942
## Test set     -13.775 38.44724 34.82500 -3.9698659 8.283390 2.435315
##                     ACF1 Theil's U
## Training set -0.10915105        NA
## Test set     -0.06905715  0.801254
accuracy(beerfit2, beer3)
##                       ME     RMSE      MAE         MPE     MAPE     MASE
## Training set   0.4761905 65.31511 54.73016  -0.9162496 12.16415 3.827284
## Test set     -51.4000000 62.69290 57.40000 -12.9549160 14.18442 4.013986
##                     ACF1 Theil's U
## Training set -0.24098292        NA
## Test set     -0.06905715  1.254009
accuracy(beerfit3, beer3)
##                     ME     RMSE  MAE        MPE     MAPE      MASE
## Training set -2.133333 16.78193 14.3 -0.5537713 3.313685 1.0000000
## Test set      5.200000 14.31084 13.4  1.1475536 3.168503 0.9370629
##                    ACF1 Theil's U
## Training set -0.2876333        NA
## Test set      0.1318407  0.298728

Base on the accuracy result, we can see the thir mothod, Seaonal naive method is the most accuracy.

non-seasonal dataset

# three different method

googfc1 <- meanf(goog200, h = 40)
googfc2 <- rwf(goog200, h = 40)
googfc3 <- rwf(goog200, h = 40, drift = T)


autoplot(subset(goog, end = 240)) +
  autolayer(googfc1, PI = F, series = "Mean") +
  autolayer(googfc2, PI = F, series = "Naive") +
  autolayer(googfc3, PI = F, series = "Drift") +
  labs(x = "day", y = "Close price (US$)",
       title = "Google Stock price forecast")

Find accuracy

googtest <- window(goog, start = 201, end = 240)

accuracy(googfc1, googtest)
##                         ME      RMSE       MAE        MPE     MAPE
## Training set -4.296286e-15  36.91961  26.86941 -0.6596884  5.95376
## Test set      1.132697e+02 114.21375 113.26971 20.3222979 20.32230
##                   MASE      ACF1 Theil's U
## Training set  7.182995 0.9668981        NA
## Test set     30.280376 0.8104340  13.92142
accuracy(googfc2, googtest)
##                      ME      RMSE       MAE       MPE      MAPE     MASE
## Training set  0.6967249  6.208148  3.740697 0.1426616 0.8437137 1.000000
## Test set     24.3677328 28.434837 24.593517 4.3171356 4.3599811 6.574582
##                     ACF1 Theil's U
## Training set -0.06038617        NA
## Test set      0.81043397  3.451903
accuracy(googfc3, googtest)
##                         ME      RMSE       MAE         MPE      MAPE
## Training set -5.998536e-15  6.168928  3.824406 -0.01570676 0.8630093
## Test set      1.008487e+01 14.077291 11.667241  1.77566103 2.0700918
##                  MASE        ACF1 Theil's U
## Training set 1.022378 -0.06038617        NA
## Test set     3.119002  0.64732736  1.709275

The best method, which have the lowest value in all summarize error term method, is drift method.

Time series cross validation

e <- tsCV(goog200, rwf, drift = T, h = 1)

# RMSE: Root mean squared error
sqrt(mean(e^2, na.rm = T))
## [1] 6.233245
sqrt(mean(residuals(rwf(goog200, drift = 1))^2, na.rm = T))
## [1] 6.168928

The RMSE from the residuals is smaller, as the corrsponding “forecast” are based on a model fitted to the entire data set, rather than being true forecasts.

A good way to choose the best forecasting model is to find the model with smallest RMSE computed using time series cross-validation.

Pipe operator

error <- goog200 %>%
  tsCV(forecastfunction = rwf, drift = T, h = 1)

error^2 %>%
  mean(na.rm = T) %>%
  sqrt()
## [1] 6.233245
error <- goog200 %>%
  rwf(drift = T) %>%
  residuals()

error^2 %>%
  mean(na.rm = T) %>%
  sqrt()
## [1] 6.168928

Example: using tsCV()

error <- tsCV(goog200, forecastfunction = naive, h = 8)

# compute the MSE value and remove missing values

mse <- colMeans(error^2, na.rm = T)

data.frame(h = 1:8, MSE = mse) %>%
  ggplot(aes(x = h, y = MSE)) +
  geom_point()

# rmse

e <- tsCV(goog200, rwf, drift = T, h = 8)

rmse <- sqrt(colMeans(error^2, na.rm = T))
rmse
##       h=1       h=2       h=3       h=4       h=5       h=6       h=7 
##  6.208148  8.578760 10.730162 12.845496 14.655090 16.082789 17.510978 
##       h=8 
## 19.150658
data.frame(h = 1:8, RESE = rmse) %>%
  ggplot(aes(x = h, y = RESE)) +
  geom_point()

3.5 Predication intervals

\(\hat{y}_{T+h|T}\pm{c}\hat{\sigma}_h\)

One step prediction intervals

# naive method residuals

residuals <- residuals(naive(goog200))

sd(residuals, na.rm = T)
## [1] 6.184487

Multiple Predication intervals

Common feature of predication intervals is that they increase in length as the forecast horizon increases. The future ahead we forecast, the more uncertainity is associated with the forecast, and thus the wider the predication intervals.

Benchmark methods

Four benchmark methods

  1. Mean forecasts: \(\hat{\sigma}_h = \hat{\sigma}\sqrt{1+1/T}\)
  2. Naive forecast: \(\hat{\sigma}_h = \hat{\sigma}\sqrt{h}\)
  3. Seasonal Naive forecast \(\hat{\sigma}_h = \hat{\sigma}\sqrt{k + 1}\) where k is the integer part of \((h-1)/m\)
  4. Drift forecasts: \(\hat{\sigma}_h = \hat{\sigma}\sqrt{h(1+h/T)}\)

Note that when h = 1 and T is large, means you need to only predict one value and you have large dataset like 200, these all give the same approximate value \(\hat{\sigma}\)

naive(goog200, h = 10)
##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 201       531.4783 523.5222 539.4343 519.3105 543.6460
## 202       531.4783 520.2267 542.7298 514.2705 548.6861
## 203       531.4783 517.6980 545.2586 510.4031 552.5534
## 204       531.4783 515.5661 547.3904 507.1428 555.8138
## 205       531.4783 513.6880 549.2686 504.2704 558.6862
## 206       531.4783 511.9900 550.9666 501.6735 561.2830
## 207       531.4783 510.4285 552.5280 499.2854 563.6711
## 208       531.4783 508.9751 553.9814 497.0627 565.8939
## 209       531.4783 507.6101 555.3465 494.9750 567.9815
## 210       531.4783 506.3190 556.6375 493.0005 569.9561
goog200[200] + sd(residuals(naive(goog200)), na.rm = T) * qnorm(0.1, lower.tail = F)
## [1] 539.404
autoplot(naive(goog200))

### Predication intervals from bootstrpped residuals

When a normal distribution for the forecast errors is an unreasonable assumption, one alternative is to use bootstrapping, which only assumes that the forecast errors are uncorrelated

Forecast error: \(e_t = y_t - \hat{y}_{t|t-1}\) which can re-rwite this as \(y_t = \hat{y}_{t|t-1} + e_t\)

==================Review

naive(goog200, bootstrap = T)
##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 201       531.4783 526.0076 537.0061 523.2307 541.0597
## 202       531.4783 523.1504 539.8016 519.8481 546.2497
## 203       531.4783 521.0374 541.6380 517.0867 550.3955
## 204       531.4783 519.1055 543.4188 514.4375 554.0778
## 205       531.4783 517.6330 544.9533 511.8730 563.3855
## 206       531.4783 516.3541 546.7031 509.7720 581.3736
## 207       531.4783 514.8630 547.9053 507.8135 581.9491
## 208       531.4783 513.7869 549.3288 505.9135 586.0329
## 209       531.4783 512.4148 550.4386 503.6003 586.3873
## 210       531.4783 510.9900 551.6743 502.1216 588.5228

In this case, they are similary (but not identical) to the predication intervals based on teh normal distribution.

Predication intervals with transformations

If a transfromation has been used, then the predication interval should be computed on the transformed scale, and the end points back-transformation to give a predication interval on teh original scale.

3.6 The forecast pacakge in R

forecast(ausbeer, h = 4)
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2010 Q3       404.6025 385.8695 423.3355 375.9529 433.2521
## 2010 Q4       480.3982 457.5283 503.2682 445.4216 515.3748
## 2011 Q1       417.0367 396.5112 437.5622 385.6456 448.4277
## 2011 Q2       383.0996 363.5063 402.6930 353.1341 413.0651

3.7 Exercises

# US net electricity generation

autoplot(usnetelec)

lambda <- BoxCox.lambda(usnetelec)
autoplot(BoxCox(usnetelec, lambda))

kpss.test(BoxCox(usnetelec, lambda))
## Warning in kpss.test(BoxCox(usnetelec, lambda)): p-value smaller than
## printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  BoxCox(usnetelec, lambda)
## KPSS Level = 1.4583, Truncation lag parameter = 3, p-value = 0.01
# Quarterly US GDP

autoplot(usgdp)

lambda <- BoxCox.lambda(usgdp)
autoplot(BoxCox(usgdp, lambda))

kpss.test(BoxCox(usgdp, lambda))
## Warning in kpss.test(BoxCox(usgdp, lambda)): p-value smaller than printed
## p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  BoxCox(usgdp, lambda)
## KPSS Level = 4.8114, Truncation lag parameter = 4, p-value = 0.01
# monthly copper prices

autoplot(mcopper)

lambda <- BoxCox.lambda(mcopper)
autoplot(BoxCox(mcopper, lambda))

kpss.test(BoxCox(mcopper, lambda))
## Warning in kpss.test(BoxCox(mcopper, lambda)): p-value smaller than printed
## p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  BoxCox(mcopper, lambda)
## KPSS Level = 6.2659, Truncation lag parameter = 6, p-value = 0.01
# monthly US domestic enplanements

autoplot(enplanements)

lambda <- BoxCox.lambda(enplanements)
autoplot(BoxCox(enplanements, lambda))

autoplot(cangas)

lambda <- BoxCox.lambda(cangas)
autoplot(BoxCox(cangas, lambda))

kpss.test(BoxCox(cangas, lambda))
## Warning in kpss.test(BoxCox(cangas, lambda)): p-value smaller than printed
## p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  BoxCox(cangas, lambda)
## KPSS Level = 7.2459, Truncation lag parameter = 6, p-value = 0.01

Because the variance is not simple increase or decrease with series, it jump during 70-90 and decrease afterward.

retail_data <- readxl::read_excel("/Users/yifeiliu/Documents/R/data/book_exercise/forecasting/retail.xlsx", skip = 1)

myts <- ts(retail_data[, "A3349873A"],
           frequency = 12, start = c(1982, 4))

autoplot(myts)

lambda <- BoxCox.lambda(myts)
autoplot(BoxCox(myts, lambda))

kpss.test(BoxCox(myts, lambda))
## Warning in kpss.test(BoxCox(myts, lambda)): p-value smaller than printed p-
## value
## 
##  KPSS Test for Level Stationarity
## 
## data:  BoxCox(myts, lambda)
## KPSS Level = 5.8234, Truncation lag parameter = 5, p-value = 0.01
# unemployment benefits in Australia
autoplot(dole)

dole_cal <- dole / monthdays(dole)
autoplot(dole_cal)

# transformation adjust for calender, population, inflation seems appropriate


# monthly accident deaths in USA
autoplot(usdeaths)

usdeaths_cal <- cbind(usdeath = usdeaths,
                       calender = usdeaths / monthdays(usdeaths))


autoplot(usdeaths_cal, facets = T)

# calender transformation seem reasonable, if we have further dataset such as ppulation. maybe population tranform is also reasonable

# quarterly clay brick production

autoplot(bricksq)

lambda <- BoxCox.lambda(bricksq)

autoplot(BoxCox(bricksq, lambda))

# all the transfromation and adjustment I know don't seem appropriate to apply on this dataset. 
beer <- window(ausbeer, start = 1992)

fc <- snaive(beer)

autoplot(fc)

res <- residuals(fc)

autoplot(res)

checkresiduals(res)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

t.test(res)
## 
##  One Sample t-test
## 
## data:  res
## t = -0.81286, df = 69, p-value = 0.4191
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -5.428075  2.285218
## sample estimates:
## mean of x 
## -1.571429

To test if the residuals are white noise we can do Portmanteau test for autocorrelation

# lag = h, h = 2*m, m = 4. fitdf = k, k = parameters in the model, model is snave, so k = 1

Box.test(res, lag = 2*4, fitdf = 1, type = "Lj")
## 
##  Box-Ljung test
## 
## data:  res
## X-squared = 32.269, df = 7, p-value = 3.621e-05

The p value is relative large, which mean the result are signiciant, so we cannot conclude that the residuals are distringuishable from a white noise series.

# WWWusage internet usage per minute

autoplot(WWWusage)

usage_mean <- rwf(WWWusage)
mean_res <- residuals(usage_mean)

checkresiduals(mean_res)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

usage_snaive <- snaive(WWWusage)
snaive_res <- residuals(usage_mean)

checkresiduals(snaive_res)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

usage_naive <- naive(WWWusage)
naive_res <- residuals(usage_naive)

checkresiduals(usage_naive)

## 
##  Ljung-Box test
## 
## data:  Residuals from Naive method
## Q* = 145.58, df = 10, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 10
# both naive or snaive is not approporate


# bricksq

autoplot(bricksq)

bricksq_naive <- naive(bricksq)
naive_res <- residuals(bricksq_naive)

checkresiduals(naive_res)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

bricksq_snaive <- snaive(bricksq)
snaive_res <- residuals(bricksq_snaive)

checkresiduals(snaive_res)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

# both are not appropriated
  1. Yes
  2. not really
  3. MAPE: mean absolute percentage error. is unit free. disadvantage, \(y_t\) cannot equal to 0. assume the unit of measurement has a meaningful zero which in some case not correct such as temperature forecasts. also put more pentalty on negative zero. MAPE not a best forecast accuracy method.
  4. Not really, you should start with simple one, such as naive, snave, rwf with drift, and use those methods as benchmark to evaluate future model.
  5. No. you need to use the model to test on out-of-sample dataset to see its accuracy level. Because fit on test set is not forecast, it’s just fit the dataset.
myts.train <- window(myts, end=c(2010,12))
myts.test <- window(myts, start=2011)


autoplot(myts) +
  autolayer(myts.train, series="Training") +
  autolayer(myts.test, series="Test")

fc <- snaive(myts.train)

accuracy(fc,myts.test)
##                     ME     RMSE      MAE       MPE      MAPE     MASE
## Training set  7.772973 20.24576 15.95676  4.702754  8.109777 1.000000
## Test set     55.300000 71.44309 55.78333 14.900996 15.082019 3.495907
##                   ACF1 Theil's U
## Training set 0.7385090        NA
## Test set     0.5315239  1.297866
checkresiduals(fc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 624.45, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24
kurtosis(residuals(fc), na.rm = T)
## [1] 1.311521
# we can see the residuals apprea to be autocorrelated and distribution also have excess kurtosis. 

# use cross validation to test accuracy

e <- tsCV(myts, snaive, h = 4)

sqrt(mean(e^2, na.rm = T))
## [1] 25.50797
resids <- myts %>%
  snaive() %>%
  residuals()

resids^2 %>%
  mean(na.rm = T) %>%
  sqrt()
## [1] 25.68062
# test sensitive to acuracy measurement

percent(length(myts.test) / length(myts))
## [1] "9.45%"

test sensitive to acuracy measurement, we can see the tes data set we only use 9.45% of original dataset. Convention of train/test should be .8/.2. So we can these how snaive model perform in these split case

myts.train <- window(myts, end = c(2007, 6))
myts.test <- window(myts, start = c(2007, 7))

autoplot(myts) +
  autolayer(myts.train, series="Training") +
  autolayer(myts.test, series="Test")

fc <- snaive(myts.train)

accuracy(fc,myts.test)
##                     ME     RMSE      MAE      MPE     MAPE     MASE
## Training set  8.969072 19.52551 15.28729 5.450456 8.272556 1.000000
## Test set     19.358333 24.14890 20.78333 6.124900 6.664263 1.359518
##                   ACF1 Theil's U
## Training set 0.7198404        NA
## Test set     0.6031662 0.4537988
checkresiduals(fc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 584.22, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24
kurtosis(residuals(fc), na.rm = T)
## [1] 1.982122

20/80 split produce accuracy result better than 10/90 split

# quarterly visitor nights for various region of australia

autoplot(visnights)

qldmetro <- (visnights[, "QLDMetro"])
   
qldmetro_tr1 <- window(qldmetro, end = c(2015, 4))
qldmetro_tr2 <- window(qldmetro, end = c(2014, 4))
qldmetro_tr3 <- window(qldmetro, end = c(2013, 4)) 

fc1 <- snaive(qldmetro_tr1)
fc2 <- snaive(qldmetro_tr2)
fc3 <- snaive(qldmetro_tr3)

accuracy(fc1, qldmetro)
##                      ME      RMSE       MAE        MPE     MAPE      MASE
## Training set 0.02006107 1.0462821 0.8475553 -0.2237701 7.976760 1.0000000
## Test set     0.56983879 0.9358727 0.7094002  4.6191866 6.159821 0.8369957
##                    ACF1 Theil's U
## Training set 0.06014484        NA
## Test set     0.09003153 0.4842061
accuracy(fc2, qldmetro)
##                     ME      RMSE       MAE        MPE     MAPE      MASE
## Training set 0.0161876 1.0735582 0.8809432 -0.2744747 8.284216 1.0000000
## Test set     0.3669560 0.8516202 0.6644561  2.8431247 6.085501 0.7542553
##                      ACF1 Theil's U
## Training set  0.066108793        NA
## Test set     -0.002021023 0.5102527
accuracy(fc3, qldmetro)
##                        ME     RMSE       MAE        MPE     MAPE      MASE
## Training set -0.007455407 1.074544 0.8821694 -0.5625865 8.271365 1.0000000
## Test set      0.411850940 1.035878 0.8800104  4.4185560 8.635235 0.9975526
##                     ACF1 Theil's U
## Training set  0.07317746        NA
## Test set     -0.21876890 0.8365723

we can see the result from training dataset 2 produce the lowest MAPE score.

autoplot(dowjones)

autoplot(rwf(dowjones, drift = T))

# this forecast are identical to extending the line drawn between the first and last observation

autoplot(rwf(dowjones, drift = T))